In [1]:

    

import ebossanalysis
import speclines
from scipy import stats


    

In [111]:

    

reload(ebossanalysis)
data = ebossanalysis.unify_emissionline_profile_readin()
emission_bstrap = ebossanalysis.unify_emissionline_profile_readin(bootstrap=True)
absorption = ebossanalysis.unify_absorptionline_profile_readin()
absorption_bstrap = ebossanalysis.unify_absorptionline_profile_readin(bootstrap=True)
corrected = ebossanalysis.corrected_velspace_flux_readin()
corrected_bstrap = ebossanalysis.corrected_velspace_flux_readin(bootstrap=True)
velmeasure = ebossanalysis.do_velocity_nonparametric()
velmeasure_bstrap = ebossanalysis.do_velocity_nonparametric(bootstrap=True)


    

In [119]:

    

#bootstrap errors
#emission_error = np.std(emission_bstrap['UNIFIEDFLUX'], axis=1)
#absorption_error = np.std(absorption_bstrap['UNIFIEDABSORPTION'], axis=1)
#bootstrap errors
emission_error = emission_bstrap['UNIFIEDDISP']
absorption_error = absorption_bstrap['UNIFIEDABSORPTION_DISP']
vpercent_lines_error = np.std(velmeasure_bstrap['FLUX_PERCENT'], axis=2)
vpercent_error = np.std(velmeasure_bstrap['UNIFIEDPERCENT'], axis=1)


    

In [225]:

    

fig, axes = subplots(figsize=(12,12), ncols=1, nrows=3)
axesall = axes.ravel()
fig.subplots_adjust(hspace=0, top=0.90, bottom=0.10)
xlimits = [-800,600]
ylimits = [0.9, 1.45]
ylimits1 = [0.9, 1.35]

#use_indices = data['INDEX']
use_indices = [0,1,2,3]
lines = data['LINES']
vel = data['VEL']
line_names = ['2366', '2396', '2613', '2626']
CM = get_cmap('jet')
line_colors = [CM(fcolor) for fcolor in (arange(4)/8.+0.00)]

for i, ax in enumerate(axesall):
    ax.tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
    ax.tick_params(axis='x', which='minor', length=4, width=2, labelsize=20)
    ax.tick_params(axis='y', which='major', length=6, width=2, labelsize=18)
    ax.tick_params(axis='y', which='minor', length=3, width=2, labelsize=18)
    ax.set_xlim(xlimits) 
    ax.plot(xlimits, [1,1], ':', color='gray', lw=3)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(2.2)

# Top panel:
ax = axesall[0]
for i,index in enumerate(use_indices):
    ax.plot(vel, data['FLUX'][index, :], color=line_colors[index], lw=4)
    ax.text(-630, 1.35-i*0.05, line_names[index], color=line_colors[index], fontsize=24)
ax.text(-750, 1.35, 'Fe II', color=line_colors[0], fontsize=22)
ax.xaxis.tick_top()
ax.set_xlabel('Velocity (km/s)', fontsize=22)
#ax.set_ylabel('Flux (normalized)', fontsize=22)
ax.text(200, 1.35, 'Observed Profile', fontsize=20)
ax.set_ylim(ylimits)
    
# Middle panel
ax = axesall[1]
for index in use_indices:
    ax.plot(vel, data['NORMFLUX'][index, :], color=line_colors[index], lw=4)
#ax.plot(vel, data['UNIFIEDFLUX'], color='black', lw=4)
setp(ax.get_xticklabels(), visible=False)
ax.set_ylabel(r'Normalized $f(\lambda)$', fontsize=22)
ax.text(200, 1.25, 'Normalized Profile', fontsize=20)
ax.set_ylim(ylimits1)
    
# Bottom panel
ax = axesall[2]
#ax.fill_between(vel, data['UNIFIEDFLUX']-emission_error*2., data['UNIFIEDFLUX']+emission_error*2., color='SkyBlue')
#ax.fill_between(vel, data['UNIFIEDFLUX']-0.02, data['UNIFIEDFLUX']+0.02, color='RoyalBlue')
ax.fill_between(vel, data['UNIFIEDFLUX']-emission_error, data['UNIFIEDFLUX']+emission_error, color='Blue', alpha=0.4)
ax.plot(vel, data['UNIFIEDFLUX'], color='Blue', lw=4)
ax.text(200, 1.25, 'Unified Profile', fontsize=20)
ax.set_xlabel('Velocity (km/s)', fontsize=22)
#ax.set_ylabel('Flux (normalized)', fontsize=22)
ax.set_ylim(ylimits1)


    

    
Out[225]:





(0.9, 1.35)






    

In [157]:

    

noffset = 100
fig, axes = subplots(figsize=(12,8), ncols=1, nrows=2)
axesall = axes.ravel()
fig.subplots_adjust(hspace=0, top=0.90, bottom=0.10)
xlimits = [-800,600]
ylimits = [0.9, 1.45]
ylimits1 = [0.9, 1.35]
ylimits2 = [0.95, 1.11]

for i, ax in enumerate(axesall):
    ax.tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
    ax.tick_params(axis='x', which='minor', length=4, width=2, labelsize=20)
    ax.tick_params(axis='y', which='major', length=6, width=2, labelsize=18)
    ax.tick_params(axis='y', which='minor', length=3, width=2, labelsize=18)
    ax.set_xlim(xlimits) 
    ax.plot(xlimits, [1,1], ':', color='gray', lw=2)
    if i == 0: ax.plot([0,0],ylimits1,':', color='black', lw=2)
    if i == 1: ax.plot([0,0],ylimits2,':', color='black', lw=2)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(2.2)
    
# Top panel
ax = axesall[0]
#ax.fill_between(vel, data['UNIFIEDFLUX']-emission_error*2., data['UNIFIEDFLUX']+emission_error*2., color='SkyBlue')
ax.fill_between(vel, data['UNIFIEDFLUX']-emission_error, data['UNIFIEDFLUX']+emission_error, color='RoyalBlue', alpha=0.7)
ax.plot(vel, data['UNIFIEDFLUX'], color='Blue', lw=4)
# Blue side
#ax.fill_between(vel[:noffset+1], data['UNIFIEDFLUX'][:noffset+1]-0.04, data['UNIFIEDFLUX'][:noffset+1]+0.04, color='SkyBlue')
#ax.fill_between(vel[:noffset+1], data['UNIFIEDFLUX'][:noffset+1]-0.02, data['UNIFIEDFLUX'][:noffset+1]+0.02, color='RoyalBlue')
#ax.plot(vel[:noffset+1], data['UNIFIEDFLUX'][:noffset+1], color='DarkBlue', lw=4)
# Red side
#ax.fill_between(vel[noffset:], data['UNIFIEDFLUX'][noffset:]-0.04, data['UNIFIEDFLUX'][noffset:]+0.04, color='LightSalmon')
#ax.fill_between(vel[noffset:], data['UNIFIEDFLUX'][noffset:]-0.02, data['UNIFIEDFLUX'][noffset:]+0.02, color='LightCoral')
#ax.plot(vel[noffset:], data['UNIFIEDFLUX'][noffset:], color='Maroon', lw=4)
# Red side flipped
#ax.fill_between(vel[noffset-20:noffset+1], data['UNIFIEDFLUX'][noffset+20:noffset-1:-1]-0.04, 
#                data['UNIFIEDFLUX'][noffset+20:noffset-1:-1]+0.04, color='LightSalmon')
#ax.fill_between(vel[noffset-20:noffset+1], data['UNIFIEDFLUX'][noffset+20:noffset-1:-1]-0.02, 
#                data['UNIFIEDFLUX'][noffset+20:noffset-1:-1]+0.02, color='LightCoral')
ax.plot(vel[noffset-20:noffset+1], data['UNIFIEDFLUX'][noffset+20:noffset-1:-1], '--', color='Red', lw=4)
ax.plot(vel[noffset:], data['UNIFIEDFLUX'][noffset:], '--', color='Red', lw=4)

ax.text(-750, 1.28, 'Unified Profile', color='Blue', fontsize=20)
ax.text(-750, 1.22, 'Symmetric Profile', color='red', fontsize=20)
ax.text(-750, 1.18, r'assuming $f(v<0)=f(|v|)$', color='red', fontsize=20)
ax.set_xlabel('Velocity (km/s)', fontsize=22)
ax.xaxis.tick_top()
ax.set_ylim(ylimits1)

ax = axesall[1]
xx = vel[noffset-20:noffset+1]
yy = data['UNIFIEDFLUX'][noffset-20:noffset+1]-data['UNIFIEDFLUX'][noffset+20:noffset-1:-1]+1.
xx1 = vel[noffset:noffset+20]
yy1 = ones(xx1.shape)
ax.plot(xx, yy, color='Blue', lw=6)
ax.fill_between(xx, yy, ones(yy.size), color='DeepSkyBlue')
ax.plot(xx1, yy1, color='Blue', lw=6)
yerr = median(emission_error[noffset-6:noffset+7])*sqrt(2)
ax.errorbar([420],[1.07],yerr=yerr,fmt='o',ms=10,lw=3,capsize=5)

#ax.fill_between(xx, yy-0.05, yresidual+0.05, color='SkyBlue')
#ax.fill_between(xx, yy-0.025, yresidual+0.025, color='RoyalBlue')
ax.text(-750, 1.08, 'Emission Excess', color='Blue', fontsize=20)
ax.text(-750, 1.065, r'at $v<0$', color='Blue', fontsize=20)
ax.set_xlabel('Velocity (km/s)', fontsize=22)
#ax.set_ylabel('Flux (normalized)', fontsize=22)
ax.set_ylim(ylimits2)
ax.text(-950, 1.06, r'Normalized $f(\lambda)$', rotation='vertical', va='bottom', fontsize=22)


    

    
Out[157]:





<matplotlib.text.Text at 0x12b626f90>






    

In [100]:

    

vel[noffset-20:noffset+1].shape, data['UNIFIEDFLUX'][noffset+20:noffset-1:-1].shape


    

    
Out[100]:





((21,), (21,))

In [161]:

    

print(speclines.FeII2366)


    

    




 Name = FeII* 2366
 Vacuum Wavelength = 2365.55 Ang
 Einstein A = 59000000.0 
 Einstein B (u->l) = 0.0 
 Einstein B (l->u) = 0.0 
 Strength = 3.08 
 Oscillator Strength = 0.0495
 Reference 1 = NIST 
 Reference 2 = NIST 

In [288]:

    

f2626_m = 2.39E-1*2600.**2#*3.52E7#/(3.52E7+2.35E8)
f2612_m = 7.17E-2*2586.**2#*1.20E8#/(8.94E7+1.20E8+6.29E7)
f2396_m = 3.59E-2*2374.**2#*2.59E8/(4.25E7+2.59E8)
f2366_m = 1.14E-1*2344.**2#*5.90E7#/(1.73E8+5.90E7+3.10E7)
f2626_s = 2.39E-1*2600.**2*3.52E7/(3.52E7+2.35E8)
f2612_s = 7.17E-2*2586.**2*1.20E8/(8.94E7+1.20E8)#+6.29E7)
f2396_s = 3.59E-2*2374.**2*2.59E8/(4.25E7+2.59E8)
f2366_s = 1.14E-1*2344.**2*5.90E7/(1.73E8+5.90E7)#+3.10E7)


    

In [289]:

    

print(f2366_s/f2626_s, f2396_s/f2626_s, f2612_s/f2626_s)
print(f2366_m/f2626_m, f2396_m/f2626_m, f2612_m/f2626_m)


    

    




(0.7568018361833609, 0.8257838969474178, 1.3055059327657133)
(0.38768184991706073, 0.12523083632492388, 0.29677792899408284)

In [301]:

    

data['FNORM']/data['FNORM'][1]


    

    
Out[301]:





array([ 0.68766044,  1.        ,  0.74631255,  1.4503818 ])

In [188]:

    

print(f2366_s/f2366_m, f2396_s/f2396_m, f2612_s/f2612_m, f2626_s/f2626_m)


    

    




(0.22433460076045625, 0.8590381426202321, 0.4406904149834741, 0.1302738712065137)

In [190]:

    

print(f2366_s/f2396_s, f2396_s/f2396_s, f2612_s/f2396_s, f2626_s/f2396_s)
print(f2366_m/f2396_m, f2396_m/f2396_m, f2612_m/f2396_m, f2626_m/f2396_m)


    

    




(0.8084403915858491, 1.0, 1.2157421556464412, 1.210970574379795)
(3.095737929204446, 1.0, 2.3698470576692707, 7.985253707045447)

In [410]:

    

data['FNORM']/data['FNORM'][3]


    

    
Out[410]:





array([ 0.47412374,  0.68947362,  0.51456282,  1.        ])

In [411]:

    

figsize(8,7)
#indices = [2,0,3]
indices = [1,0,3]
#names = array(['Fe II 2613', '2366', '2626'])
#y0 = data['FNORM'][indices]/data['FNORM'][1]
#f_s = array([f2612_s/f2396_s, f2366_s/f2396_s, f2626_s/f2396_s])
#f_m = array([f2612_m/f2396_m, f2366_m/f2396_m, f2626_m/f2396_m])
#y_s = f_s/y0
#y_m = f_m/y0
y0 = data['FNORM']/data['FNORM'][2]
#f_s = array([f2366_s/f2396_s, f2396_s/f2396_s, f2612_s/f2396_s, f2626_s/f2396_s])
#f_m = array([f2366_m/f2396_m, f2396_m/f2396_m, f2612_m/f2396_m, f2626_m/f2396_m])
f_s = array([f2366_s/f2612_s, f2396_s/f2612_s, f2612_s/f2612_s, f2626_s/f2612_s])
f_m = array([f2366_m/f2612_m, f2396_m/f2612_m, f2612_m/f2612_m, f2626_m/f2612_m])
y_s = f_s/y0
y_m = f_m/y0
x = arange(3)+1
plot([0.5,3.5],[1,1],color='Black',lw=9, alpha=0.7)


    #ax.text(-600, 1.35-i*0.05, line_names[index], color=line_colors[index], fontsize=22)
    
#plot(x, y_s, u'bs', ms=20)
#errorbar(x,y_s,yerr=0.85, color='blue', mec='blue', fmt='s', capsize=10, capthick=3, lw=3)
#plot(x, y_m, u'go', ms=20)
#errorbar(x,y_m,yerr=0.85, color='green', mec='green', fmt='o', capsize=10, capthick=3, lw=3)
xlim(0.5, 3.5)
ylim(-1,9)
for i, index in enumerate(indices):
    plot(x[i], y_s[index], u's', color=line_colors[index], ms=20)
    errorbar(x[i], y_s[index], yerr=0.85, color=line_colors[index], mec=line_colors[index], capsize=10, capthick=3, lw=3)
    plot(x[i], y_m[index], u'o', color=line_colors[index], ms=20)
    errorbar(x[i], y_m[index], yerr=0.85, color=line_colors[index], mec=line_colors[index], capsize=10, capthick=3, lw=3)

text(0.7, 7.85, r'Circle: Multiple Scattering ($\tau\gg1$)', color='black', fontsize=22)
text(0.7, 7.2, r'Square: Single Scattering ($\tau\ll1$)', color='black', fontsize=22)
ylabel(r'($\frac{W_{\rm FeII}^{\lambda}}{W_{\rm FeII}^{\lambda2396}})^{\rm e}/\frac{W_{\rm FeII}^{\lambda}}{W_{\rm FeII}^{\lambda2396}})^{\rm o}$', fontsize=22)
xticks(x, names, rotation=0, fontsize=22)
plot(x[0], y_s[2], color=line_colors[2], ms=20)


    

    
Out[411]:





[<matplotlib.lines.Line2D at 0x13366eb10>]






    

In [315]:

    

x[0], y_s[2], line_colors[2]


    

    
Out[315]:





(1, 2.11832063305962, 'red')

In [378]:

    

t2626 = 2.39E-1*2600.#*3.52E7#/(3.52E7+2.35E8)
t2612 = 7.17E-2*2586.#*1.20E8#/(8.94E7+1.20E8+6.29E7)
t2396 = 3.59E-2*2374.#*2.59E8/(4.25E7+2.59E8)
t2366 = 1.14E-1*2344.#*5.90E7#/(1.73E8+5.90E7+3.10E7)
pf2626 = 3.52E7/(3.52E7+2.35E8)
pf2612 = 1.20E8/(8.94E7+1.20E8+6.29E7)
pf2396 = 2.59E8/(4.25E7+2.59E8)
pf2366 = 5.90E7/(1.73E8+5.90E7+3.10E7)
pr2626 = 2.35E8/(3.52E7+2.35E8)
pr2612 = 8.94E7/(8.94E7+1.20E8+6.29E7)
pr2396 = 4.25E7/(4.25E7+2.59E8)
pr2366 = 1.73E8/(1.73E8+5.90E7+3.10E7)


    

In [384]:

    

tt2396 = logspace(-3,3,1000)
tt2626 = tt2396/t2396*t2626
tt2612 = tt2396/t2396*t2612
tt2366 = tt2396/t2396*t2366

beta2396 = (1.-exp(-tt2396))/tt2396
beta2626 = (1.-exp(-tt2626))/tt2626
beta2612 = (1.-exp(-tt2612))/tt2612
beta2366 = (1.-exp(-tt2366))/tt2366

PF_2396 = pf2396/(1.-pr2396*(1.-beta2396))
PF_2626 = pf2626/(1.-pr2626*(1.-beta2626))
PF_2612 = pf2612/(1.-pr2612*(1.-beta2612))
PF_2366 = pf2366/(1.-pr2366*(1.-beta2366))


    

In [402]:

    

figsize(8,7)
plot(tt2396, PF_2396, 'black', lw=3)
xscale('log')
xlim(5E-3, 1E2)
plot(tt2396, PF_2612, 'green', lw=3)
plot(tt2396, PF_2366, 'olive', lw=3)
plot(tt2396, PF_2626, 'magenta', lw=3)


    

    
Out[402]:





[<matplotlib.lines.Line2D at 0x132387fd0>]






    

In [409]:

    

figsize(8,7)
#plot(tt2396, PF_2396/PF_2396, 'black', lw=3)
xlim(5E-3, 1E2)
xscale('log')

plot(tt2396, PF_2612*t2612*2600./(PF_2396*t2396*2374.), 'olive', lw=3)
#plot([1E-2,1E2], [y0[2],y0[2]], color='Olive', lw=10, alpha=0.7)
plot(2.5E-2, y0[2], u's', color='olive', ms=15)
errorbar(2.5E-2, y0[2], yerr=0.5, color='olive', mec='olive', capsize=10, capthick=3, lw=3)
text(3.1E-2, 6.5, '2613', color='olive', fontsize=22)

plot(tt2396, PF_2366*t2366*2344./(PF_2396*t2396*2374.), 'blue', lw=3)
#fill_between([1E-2, 1E2], [y0[0]-0.3, y0[0]-0.3], [y0[0]+0.3, y0[0]+0.3], color='Blue', alpha=0.2)
#plot([1E-2,1E2], [y0[0],y0[0]], color='Blue', lw=6, alpha=0.7)
plot(4E-2, y0[0], u's', color='blue', ms=15)
errorbar(4E-2, y0[0], yerr=0.5, color='blue', mec='blue', capsize=10, capthick=3, lw=3)
text(3.1E-2, 7, '2366', color='blue', fontsize=22)

plot(tt2396, PF_2626*t2626*2600./(PF_2396*t2396*2374.), 'magenta', lw=3)
#fill_between([1E-2, 1E2], [y0[3]-0.3, y0[3]-0.3], [y0[3]+0.3, y0[3]+0.3], color='Magenta', alpha=0.2)
#plot([1E-2,1E2], [y0[3],y0[3]], color='Magenta', lw=6, alpha=0.7)
plot(5.4E-2, y0[3], u's', color='magenta', ms=15)
errorbar(5.4E-2, y0[3], yerr=0.5, color='magenta', mec='magenta', capsize=10, capthick=3, lw=3)
text(3.1E-2, 6, '2626', color='magenta', fontsize=22)

text(1E-2, 7, 'Fe II', color='blue', fontsize=22)
#text(1E-2, 5.5, 'Square: Observed Line Ratios', color='blue', fontsize=22)

xlabel(r'Optical Depth $\tau_{\rm Fe II}^{\lambda2396}$', fontsize=22)
#ylabel(r'$W^{\lambda2396}$', fontsize=22)
ylabel(r'Line Ratio $\frac{W_{\rm FeII}^{\lambda}}{W_{\rm FeII}^{\lambda2396}}$', fontsize=22)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)


    

In [ ]:

    




    

In [336]:

    

y0


    

    
Out[336]:





array([ 0.68766044,  1.        ,  0.74631255,  1.4503818 ])

In [354]:

    

1.20E8/(1.20E8+8.94E7+6.29E7), 1.20E8/(1.20E8+8.94E7)


    

    
Out[354]:





(0.4406904149834741, 0.5730659025787965)

In [383]:

    

6.29E7/1.20E8


    

    
Out[383]:





0.5241666666666667

In [404]:

    

data.dtype


    

    
Out[404]:





dtype([('LINES', '>f8', (4,)), ('VEL', '>f8', (200,)), ('FLUX', '>f8', (4, 200)), ('NORMFLUX', '>f8', (4, 200)), ('FNORM', '>f8', (4,)), ('INDEX', '>i4', (3,)), ('UNIFIEDFLUX', '>f8', (200,))])

In [427]:

    

absorption['LINES']


    

    
Out[427]:





array([ 2382.764,  2796.352,  2803.531,  2852.964,  2344.213,  2374.46 ,
        2586.649,  2600.172])

In [189]:

    

figsize(16,6)
xlimits = [-2400, 1800]
ylimits = [-0.2, 1.25]
flux_2344 = absorption['FLUX'][2]
fabs_2344 = corrected['FLUX'][2]
vel = absorption['VEL']
plot(vel+69./2., fabs_2344, '-', color='red', lw=4, drawstyle='steps', alpha=0.7)
plot(vel+69./2., flux_2344, '-', color='blue', lw=4, drawstyle='steps', alpha=1.0)
plot([500, 850], [0.135, 0.135], '-', color='blue', lw=3)
text(1000, 0.11, 'Observed', color='blue', fontsize=20)
plot([500, 850], [0.035, 0.035], '-', color='red', lw=3, alpha=0.7)
text(1000, 0.01, 'Emission Corrected', color='red', fontsize=20)
xlim(xlimits)
ylim(ylimits)
plot(xlimits, [1,1], ':', color='gray', lw=3)
plot([0,0], ylimits, '--', lw=2, color='gray')
xlabel('Velocity (km/s)', fontsize=22)
ylabel(r'Normalized $f(\lambda)$', fontsize=22)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)
text(-2200, 0.05, r'Mg II $\lambda\lambda$2796,2804', color='black', fontsize=22)


    

    
Out[189]:





<matplotlib.text.Text at 0x134ac79d0>






    

In [248]:

    

fig, axes = subplots(figsize=(16,15), ncols=1, nrows=3)
axesall = axes.ravel()
fig.subplots_adjust(hspace=0, top=0.90, bottom=0.10)
xlimits = [-2400, 1800]
ylimits = [-0.2, 1.25]

indices = [0,7,2]
vel = absorption['VEL']
for i, index in enumerate(indices):
    ax = axesall[i]
    this_flux = absorption['FLUX'][index]
    this_fabs = corrected['FLUX'][index]
    ax.plot(vel+69./2., this_fabs, '-', color='red', lw=4, drawstyle='steps', alpha=0.7)
    ax.plot(vel+69./2., this_flux, '-', color='blue', lw=4, drawstyle='steps', alpha=1.0)
    if (i == 1):
        ax.plot([500, 850], [0.135, 0.135], '-', color='blue', lw=3)
        ax.text(1000, 0.11, 'Observed', color='blue', fontsize=20)
        ax.plot([500, 850], [0.035, 0.035], '-', color='red', lw=3, alpha=0.7)
        ax.text(1000, 0.01, 'Emission Corrected', color='red', fontsize=20)
    ax.set_xlim(xlimits)
    ax.set_ylim(ylimits)
    ax.plot(xlimits, [1,1], ':', color='gray', lw=3)
    ax.plot([0,0], ylimits, '--', lw=2, color='gray')
    if (i == 0):
        ax.text(-2200, 0.05, r'Fe II $\lambda$2374 & $\lambda$2383', color='black', fontsize=22)
    if (i == 1):
        ax.text(-2200, 0.05, r'Fe II $\lambda$2587 & $\lambda$2600', color='black', fontsize=22)
    if (i == 2):
        ax.text(-2200, 0.05, r'Mg II $\lambda\lambda$2796,$\,$2804', color='black', fontsize=22)
        
xlabel('Velocity (km/s)', fontsize=22)
ylabel(r'Normalized $f(\lambda)$', fontsize=22)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)


    

In [246]:

    

xlimits = [-800,600]
ylimits = [0.75, 1.075]

#use_indices = data['INDEX']
#use_indices = [3,5,6,4,7,0,2,1]
use_indices = (array([5, 6, 4, 7, 0, 3, 2, 1]))[::-1]
lines = absorption['LINES']
vel = absorption['VEL']
line_names = [str(int(rint(lines[index]))) for index in arange(8)]

figsize(18,6)
# Middle panel
CM = get_cmap('jet')
line_colors = [CM(fcolor) for fcolor in arange(8)/8.+0.]
for i, index in enumerate(use_indices):
    plot(vel, 1.-(1.-absorption['FLUX'][index, :])/velmeasure['TFLUX'][index]*0.55, color=line_colors[i], lw=4)
    text(-740, 0.90-i*0.018, line_names[index], color=line_colors[i], fontsize=22)

plot(xlimits, [1,1], ':', color='gray', lw=3)
plot([0,0], ylimits, '--', lw=2, color='gray')

xlabel('Velocity (km/s)', fontsize=22)
ylabel(r'Normalized $f(\lambda)$', fontsize=22)
text(200, 0.81, 'Observed Absorption Profile', fontsize=20)
text(200, 0.79, '(Normalized)', fontsize=20)
xlim(xlimits)
ylim(ylimits)
#fill_between(vel, absorption['UNIFIEDABSORPTION']-absorption_error, 
#                absorption['UNIFIEDABSORPTION']+absorption_error, color='Red')


    

    
Out[246]:





(0.75, 1.075)






    

In [188]:

    

figsize(16,6)
#xlimits = [-1800, 3500]
xlimits = [-2400, 1800]
ylimits = [-0.2, 1.25]
flux_2344 = absorption['FLUX'][7]
fabs_2344 = corrected['FLUX'][7]
vel = absorption['VEL']
plot(vel+69./2., fabs_2344, '-', color='red', lw=4, drawstyle='steps', alpha=0.7)
plot(vel+69./2., flux_2344, '-', color='blue', lw=4, drawstyle='steps', alpha=1.0)
plot([500, 850], [0.135, 0.135], '-', color='blue', lw=3)
text(1000, 0.11, 'Observed', color='blue', fontsize=20)
plot([500, 850], [0.035, 0.035], '-', color='red', lw=3, alpha=0.7)
text(1000, 0.01, 'Emission Corrected', color='red', fontsize=20)
xlim(xlimits)
ylim(ylimits)
plot(xlimits, [1,1], ':', color='gray', lw=3)
plot([0,0], ylimits, '--', lw=2, color='gray')
xlabel('Velocity (km/s)', fontsize=22)
ylabel(r'Normalized $f(\lambda)$', fontsize=22)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)
text(-2200, 0.11, r'Fe II $\lambda$2587', color='black', fontsize=22)
text(-2200, 0.02, r'Fe II $\lambda$2600', color='black', fontsize=22)


    

    
Out[188]:





<matplotlib.text.Text at 0x134a7be90>






    

In [249]:

    

fig, axes = subplots(figsize=(12,12), ncols=1, nrows=3)
axesall = axes.ravel()
fig.subplots_adjust(hspace=0, top=0.90, bottom=0.10)
xlimits = [-800,600]
ylimits = [-0.1, 1.15]
ylimits1 = [0.75, 1.075]

#use_indices = data['INDEX']
use_indices = [3,5,6,4,7,0,2,1]
lines = absorption['LINES']
vel = absorption['VEL']
line_names = [str(int(rint(lines[index]))) for index in arange(8)]
#line_names = ['2853', '2374', '2587', '2344', '2600', '2382', '2803', '2796']
#line_names = [2796', '2803', '2382', '2600', '2344', '2587', '2374', '2853']
#line_colors = ['blue', 'green', 'olive', 'magenta']

for i, ax in enumerate(axesall):
    ax.tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
    ax.tick_params(axis='x', which='minor', length=4, width=2, labelsize=20)
    ax.tick_params(axis='y', which='major', length=6, width=2, labelsize=18)
    ax.tick_params(axis='y', which='minor', length=3, width=2, labelsize=18)
    ax.set_xlim(xlimits) 
    ax.plot(xlimits, [1,1], ':', color='gray', lw=3)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(2.2)

# Top panel:
CM = get_cmap('gist_heat')
line_colors = [CM(fcolor) for fcolor in arange(8)/8./1.5+0.1]
ax = axesall[0]
for i, index in enumerate(use_indices):
    ax.plot(vel, absorption['FABS'][index], color=line_colors[i], lw=2)
    ax.text(-740, 0.70-i*0.084, line_names[index], color=line_colors[i], fontsize=22)

    
#ax.text(-750, 1.35, 'Fe II', color=line_colors[0], fontsize=22)
ax.xaxis.tick_top()
ax.set_xlabel('Velocity (km/s)', fontsize=22)
#ax.set_ylabel('Flux (normalized)', fontsize=22)
ax.text(150, 0.2, 'Emission Corrected', fontsize=20)
ax.text(150, 0.1, 'Absorption Profile', fontsize=20)
ax.set_ylim(ylimits)
    
# Middle panel
ax = axesall[1]
for i, index in enumerate(use_indices):
    color = CM(i/8./1.5+0.2)
    ax.plot(vel, absorption['NORMFABS'][index, :], color=color, lw=3)
#ax.plot(vel, data['UNIFIEDFLUX'], color='black', lw=4)
setp(ax.get_xticklabels(), visible=False)
ax.set_ylabel(r'Normalized $f(\lambda)$', fontsize=22)
ax.text(200, 0.8, 'Normalized Profile', fontsize=20)
ax.set_ylim(ylimits1)
    
# Bottom panel
ax = axesall[2]
#ax.fill_between(vel, absorption['UNIFIEDABSORPTION']-absorption_error*2, 
#                absorption['UNIFIEDABSORPTION']+absorption_error*2, color='Bisque')
ax.fill_between(vel, absorption['UNIFIEDABSORPTION']-absorption_error, 
                absorption['UNIFIEDABSORPTION']+absorption_error, color='LightSalmon')
ax.plot(vel, absorption['UNIFIEDABSORPTION'], color='Maroon', lw=4)
ax.text(200, 0.8, 'Unified Profile', fontsize=20)
ax.set_xlabel('Velocity (km/s)', fontsize=22)
#ax.set_ylabel('Flux (normalized)', fontsize=22)
ax.set_ylim(ylimits1)


    

    
Out[249]:





(0.75, 1.075)






    

In [258]:

    

# Symmetric Gaussian profile
xtmp = linspace(-800., 600, num=100)
ytmp = 1.+stats.norm.pdf((xtmp+0)/sqrt(69.**2+50.**2))*0.6

figsize(14,6)
yabs = 1.+(1.-absorption['UNIFIEDABSORPTION'])*1.3#-0.005
#yemi = data['UNIFIEDFLUX']#+0.005
yemi = data['UNIFIEDFLUX']
plot(vel, yemi, color='blue', lw=3)#, drawstyle='steps')
fill_between(vel, yemi-emission_error, yemi+emission_error, color='RoyalBlue', alpha=0.7)
plot(vel, yabs, color='red', lw=3)#, drawstyle='steps')
fill_between(vel, yabs-absorption_error, yabs+absorption_error*1.25, color='LightSalmon',alpha=0.7)
xlim(-800,600)
ylim(0.9, 1.35)
plot(xlimits, [1,1], ':', color='black', lw=2)
xlabel('Velocity (km/s)', fontsize=22)
ylabel(r'Normalized $f(\lambda)$', fontsize=22)
plot([0,0],[0.9,1.35],':', color='black', lw=2)
plot(xtmp, ytmp, '--', lw=6, color='green')
text(-750, 1.3, 'Unified Absorption Profile (Flipped)', color='Red', fontsize=22)
text(-750, 1.27, 'Unified Emission Profile', color='Blue', fontsize=22)
text(-750, 1.24, r'Gaussian ($\sigma=85\,$km/s)', color='Green', fontsize=22)


    

    
Out[258]:





<matplotlib.text.Text at 0x13f14a850>






    

In [507]:

    

repr(int(2796.0))


    

    
Out[507]:





'2796'

In [651]:

    

sqrt(69**2+60**2)


    

    
Out[651]:





91.438503924769023

In [259]:

    

figsize(9,9)
xlimits = [-1, 8]
ylimits = [150, -490]
indices = array([5, 6, 4, 7, 0, 3, 2, 1]) #to be ordered by FR
lines = velmeasure['LINES']
names = array([str(int(rint(lines[index]))) for index in arange(8)])
v80 = velmeasure['FLUX_PERCENT'][14,:]
v80_error = vpercent_lines_error[14,:]

v50 = velmeasure['FLUX_PERCENT'][8,:]
v50_error = vpercent_lines_error[8,:]

v20 = velmeasure['FLUX_PERCENT'][2,:]
v20_error = vpercent_lines_error[2,:]

v80_corrected = velmeasure['UNIFIEDPERCENT'][14]
v80_corrected_error = vpercent_error[14]
v50_corrected = velmeasure['UNIFIEDPERCENT'][8]
v50_corrected_error = vpercent_error[8]
v20_corrected = velmeasure['UNIFIEDPERCENT'][2]
v20_corrected_error = vpercent_error[2]

xx = arange(8)
#plot(xlimits, [v80_corrected, v80_corrected], lw=2, color='Blue')
fill_between(xlimits, [v80_corrected, v80_corrected]-v80_corrected_error, [v80_corrected, v80_corrected]+v80_corrected_error,
             color='RoyalBlue', alpha=0.7)
plot(xx, v80[indices], 'o', ms=10, color='blue')
errorbar(xx, v80[indices], v80_error[indices], ms=15, color='blue', elinewidth=2, capsize=5, capthick=2)
for i, index in enumerate(indices):
    text(xx[i]-0.3, v80[index]-30., names[index], fontsize=18)
text(xx[0]-0.5, v80_corrected-160., r'$v_{80}$ (80% of absorption is at $v>v_{80}$)', fontsize=22, color='blue')
text(xx[0]-0.4, v80_corrected-280., 'Squares: Observed', fontsize=20, color='black')
text(xx[0]-0.4, v80_corrected-250., 'Horizontal bands: Emission corrected', fontsize=20, color='black')

    
#plot(xlimits, [v50_corrected, v50_corrected], lw=2, color='green', alpha=0.7)
fill_between(xlimits, [v50_corrected, v50_corrected]-v50_corrected_error, [v50_corrected, v50_corrected]+v50_corrected_error,
             color='green', alpha=0.7)
plot(xx, v50[indices], 'o', ms=10, color='green')
errorbar(xx, v50[indices], v50_error[indices], ms=15, color='green', elinewidth=2, capsize=5, capthick=2)
text(xx[0]-0.5, v50_corrected-40., r'$v_{50}$', fontsize=24, color='green')

#plot(xlimits, [v20_corrected, v20_corrected], lw=2, color='Salmon', alpha=0.7)
fill_between(xlimits, [v20_corrected, v20_corrected]-v50_corrected_error, [v20_corrected, v20_corrected]+v20_corrected_error,
             color='red', alpha=0.7)
plot(xx, v20[indices], 'o', ms=10, color='red')
errorbar(xx, v20[indices], v20_error[indices], ms=15, color='red', elinewidth=2, capsize=5, capthick=2)
text(xx[0]-0.5, v20_corrected-40., r'$v_{20}$', fontsize=24, color='red')

xlim(xlimits)
ylim(ylimits)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)
xticks(xx, names[indices], rotation=0, fontsize=20)
xlabel('Lines', fontsize=22)
ylabel('Velocity (km/s)', fontsize=22)


    

    
Out[259]:





<matplotlib.text.Text at 0x13b267d90>






    

In [845]:

    

lines = velmeasure['LINES']
indices = array([5, 6, 4, 7, 0, 3, 2, 1]) #to be ordered by FR
names = array([str(int(rint(lines[index]))) for index in arange(8)])
observed_ratio = velmeasure['TFLUX'][indices]/velmeasure['TFLUX'][5]
corrected_ratio = velmeasure['TFABS'][indices]/velmeasure['TFABS'][5]
qso_ratio = velmeasure['TFABS'][indices]/velmeasure['TFABS'][5]
xx = arange(8)
plot(xx, observed_ratio/qso_ratio, 'o', ms=15, color='blue')
errorbar(xx, observed_ratio/qso_ratio, 0.04, ms=15, color='blue')
plot(xx, corrected_ratio/qso_ratio, 's', ms=15, color='green')
errorbar(xx, corrected_ratio/qso_ratio, 0.04, ms=15, color='green')
xlimits = [-1, 8]
ylimits = [0.1, 1.1]
xlim(xlimits)
ylim(ylimits)
tick_params(axis='x', which='major', length=8, width=2, labelsize=20)
tick_params(axis='y', which='major', length=8, width=2, labelsize=20)
xticks(xx, names[indices], rotation=0, fontsize=20)
#ylabel('Velocity (km/s)', fontsize=22)
#ylabel(r'Line Ratio $\frac{W_{\rm FeII}^{\lambda}}{W_{\rm FeII}^{\lambda2374}}$', fontsize=22)
ylabel(r'($\frac{W^{\lambda}}{W^{\lambda2374}})^{\rm gal}/\frac{W^{\lambda}}{W^{\lambda2374}})^{\rm qso}$', fontsize=22)
for i, index in enumerate(indices):
    text(xx[i]-0.3, (observed_ratio/qso_ratio)[i]-0.1, names[index], fontsize=18)


    

In [62]:

    

absorption_bstrap.dtype


    

    
Out[62]:





dtype([('LINES', '>f8', (8,)), ('VEL', '>f8', (200,)), ('FLUX', '>f8', (8, 200, 500)), ('FABS', '>f8', (8, 200, 500)), ('NORMFABS', '>f8', (8, 200)), ('FNORM', '>f8', (8,)), ('INDEX', '>i4', (8,)), ('UNIFIEDABSORPTION', '>f8', (200, 500)), ('UNIFIEDEMISSION', '>f8', (200, 500)), ('COEFF', '>f8', (8, 2, 500)), ('FABS_2374', '>f8', (200, 500))])

In [63]:

    

emission_bstrap.dtype


    

    
Out[63]:





dtype([('LINES', '>f8', (4,)), ('VEL', '>f8', (200,)), ('FLUX', '>f8', (4, 200, 500)), ('NORMFLUX', '>f8', (4, 200, 500)), ('FNORM', '>f8', (4, 500)), ('INDEX', '>i4', (3,)), ('UNIFIEDFLUX', '>f8', (200, 500))])

In [78]:

    

yerr = median(absorption_error[noffset-10:noffset+10])*sqrt(2)
print yerr


    

    




0.0048251442501

In [109]:

    

corrected_bstrap.dtype


    

    
Out[109]:





dtype([('LINES', '>f8', (12,)), ('VEL', '>f8', (200,)), ('FLUX', '>f8', (12, 200))])

In [122]:

    

velmeasure_bstrap.dtype


    

    
Out[122]:





dtype([('LINES', '<f8', (8,)), ('PERCENT', '<f8', (17,)), ('FLUX_PERCENT', '<f8', (17, 8, 500)), ('FABS_PERCENT', '<f8', (17, 8, 500)), ('UNIFIEDPERCENT', '<f8', (17, 500)), ('TFLUX', '<f8', (8, 500)), ('TFABS', '<f8', (8, 500))])

In [123]:

    

velmeasure.dtype


    

    
Out[123]:





dtype([('LINES', '<f8', (8,)), ('PERCENT', '<f8', (17,)), ('FLUX_PERCENT', '<f8', (17, 8)), ('FABS_PERCENT', '<f8', (17, 8)), ('UNIFIEDPERCENT', '<f8', (17,)), ('TFLUX', '<f8', (8,)), ('TFABS', '<f8', (8,))])

In [148]:

    

median(v50_error)


    

    
Out[148]:





5.2339019076103215

In [175]:

    

absorption['LINES']


    

    
Out[175]:





array([ 2382.764,  2796.352,  2803.531,  2852.964,  2344.213,  2374.46 ,
        2586.649,  2600.172])

In [ ]: